Flight Planning

Environmental Data

In this file, the goal is to quantify the historical cloud cover and wind speeds in the planned flight boxes and to assess the implications for BioSCape operations.

To assess potential cloud cover during the BioSCape campaign, here we use cloud cover data from the MODIS aqua (https://developers.google.com/earth-engine/datasets/catalog/MODIS_061_MOD09GA) and terra (https://developers.google.com/earth-engine/datasets/catalog/MODIS_061_MOD09GA) daily surface reflectance products with a 1km resolution. The QA metadata for these products assigns each raster cell one of 4 cloud categories: Clear, Cloudy, Mixed, or “Not set, assumed clear”. Here, we consider any raster cell categorized as “cloudy” to be cloud covered and other categories to be cloud-free (Figure 1). To assess potential wind speeds, we used wind speed data from ERA5. Wind data are hourly and have a 0.25 degree resolution.

Show the code

# Load required packages

# Load packages
  library(rgee)
  library(targets)
  library(sf)
  library(terra)
  library(raster)
  library(tidyverse)
  library(lubridate)
  library(leaflet)
  library(ggplot2)
  library(ggpubr)
  library(leafem)
  library(plotly)
  library(tidyterra)




#Load required data

  # get domain
    domain <- st_read("data/output/domain.gpkg",
                      quiet = TRUE)

    domain_sf <- domain
  
  # get flight boxes
    boxes <- st_read("data/flight_planning/v2_20230718_G3_AVIRISNG_PRISM_boxes.gpkg",
                     quiet = TRUE)
    
    boxes$id <- 1:nrow(boxes) # need a unique ID to make things easier
    
    
    
    #need to set up an ordered id.  upper left to lower right

    #we'll order by centroid distance relative to the upper left bounding box of the domain

      boxes %>%
        st_centroid() -> box_centroids
      
      
      domain%>%
        st_transform(crs = st_crs(boxes))%>%
      st_bbox()%>%
      st_as_sfc() %>%
        st_cast("POINT") %>%
        st_as_sf()-> point
      
      box_dist <-st_distance(box_centroids, point[4,])
      
      boxes$ordered_id[order(box_dist)] <-1:nrow(boxes) 

      boxes_sf <- boxes
    
  # Download table from drive (to see the code underlying this or to update the data, see the file "R/mock_flights_earth_engine.R")
    
    
    cloud_table <- read.csv("data/test_cloud_stats.csv")

  # modis clouds

    cloud_table %>%
      mutate(year = year(date),
             month = month(date),
             day = day(date),
             day_of_year = yday(date)) -> cloud_table
  
  #Load era5 wind data
    
    era5_wind_table <- readRDS("data/output/era_wind_weighted.RDS")

Show the code

#Load in the correct projection (for some reason this is handled incorrectly otherwise)

  nasa_proj <- "+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +R=6371007.181 +units=m +no_defs"

#Load the example layers
    
  list.files(path = "data/flight_planning/",
             pattern = "example_cloud_cover",
             full.names = TRUE) %>%
    rast() -> cloud_examples
  
# Generate the bounding box used
    domain_plus_boxes <-
    st_union(domain_sf%>%st_transform(crs = st_crs(boxes_sf)), boxes_sf) %>%
      st_bbox() %>%
      st_as_sfc()

  crs(cloud_examples) <- nasa_proj

  
  ggplot()+
    geom_spatraster(data = cloud_examples[[1]])+
  geom_sf(data = domain)

Figure 1. Binary MODIS cloud data. Clouds are in white, non-clouds in blue.

Show the code

  cloud_examples %>%
    project(y = terra::crs(domain_plus_boxes,proj=TRUE)) %>%
    crop(y = vect(domain_plus_boxes)) -> cloud_examples

c1 <-  
  ggplot()+
    geom_spatraster(data = cloud_examples[[1]])+
    scale_fill_gradient(low = "light blue",high = "white")+
    #scale_fill_manual(values = c("light blue","white"))+
    scale_x_continuous(expand = (c(0,0)))+
    scale_y_continuous(expand = c(0,0))+
    geom_sf(data = domain, fill=NA)+
    theme(legend.position = "none")

c2 <-  
  ggplot()+
    geom_spatraster(data = cloud_examples[[2]])+
    scale_fill_gradient(low = "light blue",high = "white")+
    #scale_fill_manual(values = c("light blue","white"))+
    scale_x_continuous(expand = (c(0,0)))+
    scale_y_continuous(expand = c(0,0))+
    geom_sf(data = domain, fill=NA)+
    theme(legend.position = "none")

c3 <-  
  ggplot()+
    geom_spatraster(data = cloud_examples[[3]])+
    scale_fill_gradient(low = "light blue",high = "white")+
    #scale_fill_manual(values = c("light blue","white"))+
    scale_x_continuous(expand = (c(0,0)))+
    scale_y_continuous(expand = c(0,0))+
    geom_sf(data = domain, fill=NA)+
    theme(legend.position = "none")

ggarrange(c1,c2,c3,
          common.legend = FALSE,
          ncol = 1)

Figure 1. Binary MODIS cloud data. Clouds are in white, non-clouds in blue.

Cloud Cover and Wind Speed

To visualize spatial patterns of cloud cover and wind speed, we calculated the averages for each raster cell (Figure 2). For cloud cover, to took the mean value across days (October-December) and years (2000-present). For wind speed, we took median values. We also include median wind speed estimates from FEWSNET (https://developers.google.com/earth-engine/datasets/catalog/NASA_FLDAS_NOAH01_C_GL_M_V001), which are monthly estimates at ~11km resolution.

Code

#Pull in other bioscape layers

  flights <-
  boxes_sf %>%
  st_transform(crs = st_crs(4326)) #-> flights

  team_requests <- st_read("data/manual_downloads/BIOSCAPE_proposed/20221014_team_polygons.gpkg",quiet = TRUE) %>%
      st_transform(crs = st_crs(4326))

  domain_sf %>%
      st_transform(crs = st_crs(4326)) -> domain_wgs84
  
  mean_cloud_cover <- raster("data/output/mean_cloud_cover.tif")
  

#Make a palette
  pal <- colorNumeric(palette = colorRamp(c("white", "blue"), interpolate = "spline"),
                      domain = unique(flights$prop_clear))
  
  pal2 <- colorNumeric(palette = "Blues",
                      domain = 0:1, #unique(values(mean_cloud_cover)),
                      reverse = TRUE)

#Get era wind data  
  median_era5_wind_speed <- raster("data/output/median_era5_speed.tif")

#Get fewsnet wind data  
  median_fewsnet_wind_speed <- raster("data/output/median_wind_fewsnet.tif")

#Make a palette
  pal_era5 <- colorNumeric(palette = colorRamp(c("white", "magenta"), interpolate = "spline"),
                      domain = unique(c(values(median_era5_wind_speed),values(median_fewsnet_wind_speed))),na.color = NA)

#Make labels
labels <- sprintf(as.character(flights$prop_clear)) %>%
  lapply(htmltools::HTML)


  boxes_sf %>%
  st_transform(crs = st_crs(4326)) %>%
leaflet() %>%
  addProviderTiles("Esri.NatGeoWorldMap", group = "NatGeo") %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "World Imagery") %>%
  addMapPane("flights", zIndex = 420) %>%
  addMapPane("requests", zIndex = 410) %>% 
  addRasterImage(x = mean_cloud_cover,
                 group = "Clouds",
                 colors = pal2,
                 opacity = 1) %>%
      addRasterImage(x = median_era5_wind_speed,
                 group = "Wind E",
                 colors = pal_era5,
                 opacity = 1) %>%
      addRasterImage(x = median_fewsnet_wind_speed,
                 group = "Wind F",
                 colors = pal_era5,
                 opacity = 1) %>%
     addPolygons(stroke = TRUE,
              group = "Flights",
              color = "black",
              opacity = 1,
              weight = 1,
              label = ~as.character(ordered_id),
              labelOptions = labelOptions(noHide = T,
                                          textOnly = T,
                                          textsize = 3),
              options = pathOptions(pane = "flights"),
              fill = FALSE)%>%
      addPolygons(data = team_requests%>%
                st_zm(drop = T, what = "ZM"),
                  stroke = TRUE,
                  color = "black",
                  group = "Requests",
              options = pathOptions(pane = "requests"),
              fill = FALSE,
              weight = 1)%>%  
    addMouseCoordinates() %>%
    #addImageQuery(sampling_options_wgs84, type="mousemove", layerId = "park_name") %>%
  leaflet::addLegend(position = "bottomright",
            pal = pal2,
            values = 0:1,
            opacity = 1,
            title = "Cloud Cover") %>%
    leaflet::addLegend(position = "bottomleft",
            pal = pal_era5,
            values = min(na.omit(c(values(median_era5_wind_speed),values(median_fewsnet_wind_speed)))):max(na.omit(c(values(median_era5_wind_speed),values(median_fewsnet_wind_speed)))),
            opacity = 1,
            title = "Wind (m/s)",) %>%
    addLayersControl(
    baseGroups = c("World Imagery","NatGeo"),
    overlayGroups = c("Flights","Requests","Clouds","Wind E", "Wind F"),
    options = layersControlOptions(collapsed = FALSE),
    position = "topright") %>%
  hideGroup(c("Requests","Wind E", "Wind F"))

Figure 2. Mean cloud cover. Whiter raster cells are cloudier, bluer cells are clearer. Boxes shown are flight boxes, labeled with a unique ID. Polygons are the requested sampling regions. Wind E = Median Wind Speed from ERA5 data.Wind F = Median Wind Speed from FEWSNEWS data

Median Monthly Wind Speed and Cloud Cover

We also did spatiotemporal aggregations of mean wind speed by month and flight box (Figure 4). In other words, we aggregated both spatially (all cells within the domain) and temporally (all days within the given month across all years).

Code

# mean monthly wind speed
    era5_wind_table %>%
      group_by(ID, month)%>%
      filter(month != 9)%>%
      summarize(mean_cc = median(na.omit(wgt_median_wind_speed)))%>%
      inner_join(x = boxes_sf,
                 by = c("id"="ID"))%>%
      mutate(month = lubridate::month(month,label = TRUE)) %>%
      ggplot(mapping = aes(fill = mean_cc))+
      geom_sf()+
      geom_sf(data = domain_sf,
              inherit.aes = FALSE,fill=NA)+
      scale_fill_gradient(low = "white",high = "magenta")+
      facet_wrap(~month,ncol = 1)+
      geom_sf_text(aes(label = round(mean_cc,digits = 2)),size=2)+
      labs(fill="Mean\nWind\nSpeed\n(m/s)")+
      xlab(NULL)+
      ylab(NULL) -> monthly_wind

    cloud_table %>%
      group_by(id, month)%>%
      filter(month != 9)%>%
      summarize(mean_cc = mean(na.omit(mean)))%>%
      inner_join(x = boxes_sf,
                 by = c("id"="id"))%>%
      mutate(month = lubridate::month(month,label = TRUE)) %>%
      ggplot(mapping = aes(fill = mean_cc))+
      geom_sf()+
      geom_sf(data = domain_sf,
              inherit.aes = FALSE,fill=NA)+
      scale_fill_gradient(low = "light blue",high = "white")+
      facet_wrap(~month,ncol = 1)+
      geom_sf_text(aes(label = round(mean_cc,digits = 2)),size=2)+
      labs(fill="Mean\nCloud\nCover")+
      xlab(NULL)+
      ylab(NULL) -> monthly_cloud

  ggarrange(monthly_cloud,monthly_wind,legend = "bottom")  

Figure 2. The full domain (light grey) and planned flight boxes, colored by mean monthly wind speed. Numbers on flight boxes refer to mean wind speed.

Code

#Combine wind and cloud info

cloud_table %>%
  left_join(y = boxes_sf %>%
      st_drop_geometry())%>%
  mutate(mean_clear_skies = 1-mean)%>%
  rename(mean_cloud_cover = mean)%>%
  left_join(y = 
era5_wind_table %>%
  left_join(y = boxes_sf %>%
      st_drop_geometry(),by=c("ID"="id"))%>%
  group_by(ordered_id,day) %>%
  summarize(median_wind_speed = median(na.omit(wgt_median_wind_speed))))%>%
  pivot_longer(cols = c(mean_cloud_cover,median_wind_speed))%>%
  ungroup%>%
  mutate(ordered_id = as.factor(ordered_id))%>%
  ggplot(mapping = aes(y=value,x=ordered_id,color=name))+
  geom_violin()+facet_wrap(~name,scales = "free")

Cloud cover over time

To visualize temporal patterns in cloud cover, we calculated the mean cloud cover for each flight box (Figure 3).

Code

##| out.height: 30%


  #ID by day of year
  
    cloud_table %>%
        mutate(date = as.Date(day_of_year,origin="2022-12-31")) %>% #this line is a bit of a hack
        mutate(start_of_month = floor_date(as_date(date),unit = "month"),
               month_label = month(start_of_month,label = TRUE),
               julian_label = yday(start_of_month),
               day_of_month = mday(date),
               md_label = paste(month_label,"-",day_of_month,sep = "")) %>%
      group_by(id,day_of_year) %>%
      left_join(y = boxes_sf %>%
      st_drop_geometry())%>%
  filter(month %in% c(10,11,12)) -> temp
      

    temp %>%
      ggplot() +
      geom_tile(mapping = aes(x = day_of_year,
                              y = ordered_id,
                              fill = mean))+
      scale_fill_gradient(low = "sky blue",
                          high = "white")+
        scale_x_continuous(breaks = temp$julian_label,
                           labels = temp$month_label)+
      labs(fill="Mean\nCloud\nCover",
           y="Flight Box ID",
           x = "Day of Year")+
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              aspect.ratio = 1)+
      facet_wrap(~year)

    
    
      
    # for(i in unique(temp$year)){
    #   
    #   temp %>%
    #     filter(year == i)%>%
    #     ggplot() +
    #     geom_tile(mapping = aes(x = day_of_year,
    #                             y = ordered_id,
    #                             fill = mean))+
    #     scale_fill_gradient(low = "sky blue",
    #                         high = "white")+
    #       scale_x_continuous(breaks = temp%>%
    #                            filter(year==i) %>%
    #                            pull(julian_label),
    #                          labels = temp%>%
    #                            filter(year==i)%>%
    #                            pull(month_label),
    #                          limits = c(min(temp$julian_label),
    #                                     max(temp$julian_label)))+
    #     labs(fill="Mean\nCloud\nCover",
    #          y="Flight Box ID",
    #          x = "Day of Year",
    #          title = i) -> tt
    #   
    #   plot(tt)
    #   
    #   
    #   
    # }
    # 

Figure 3. Rows represent different flight boxes (see Fig 2.), columns different days.

Wind Speed over time

Code

##| out.height: 30%

  #ID by day of year
  
    era5_wind_table %>%
        mutate(date = as.Date(doy,origin="2022-12-31")) %>% #this line is a bit of a hack
        #mutate(date = as_date(time))%>%
        mutate(start_of_month = floor_date(date,unit = "month"),
               month_label = month(start_of_month,label = TRUE),
               julian_label = yday(start_of_month),
               day_of_month = mday(date),
               day_of_year = doy,
               md_label = paste(month_label,"-",day_of_month,sep = "")) %>%
      group_by(ID,doy,day_of_year,month,month_label,julian_label,year) %>%
      summarise(median = median(wgt_median_wind_speed))%>%
      left_join(y = boxes_sf %>%
      st_drop_geometry(),by = c("ID"="id"))%>%
      filter(month %in% c(10,11,12))-> temp

    temp %>%
      ggplot() +
      geom_tile(mapping = aes(x = day_of_year,
                              y = ordered_id,
                              fill = median))+
      scale_fill_gradient(low = "white",
                          high = "magenta")+
      facet_wrap(~year)+
        scale_x_continuous(breaks = temp$julian_label,
                           labels = temp$month_label)+
      labs(fill="Median\nWind\nSpeed\n(m/s)",
           y="Flight Box ID",
           x = "Day of Year")+
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

    # for(i in unique(temp$year)){
    # 
    # 
    #   temp %>%
    #     filter(year == i)%>%
    #     ggplot() +
    #     geom_tile(mapping = aes(x = day_of_year,
    #                             y = ordered_id,
    #                             fill = median))+
    #     scale_fill_gradient(low = "white",
    #                         high = "magenta")+
    #       scale_x_continuous(breaks = temp%>%
    #                            filter(year==i) %>%
    #                            pull(julian_label),
    #                          labels = temp%>%
    #                            filter(year==i)%>%
    #                            pull(month_label),
    #                          limits = c(min(temp$julian_label),
    #                                     max(temp$julian_label)))+
    #     labs(fill="Median\nWind\nSpeed",
    #          y="Flight Box ID",
    #          x = "Day of Year",
    #          title = i) -> tt
    #   
    #   plot(tt)
    #   
    #   
    #   
    # }

Figure 4. Rows represent different flight boxes (see Fig 2.), columns different days.

Proportion of clear days for each flight box

Here, we define a “clear” day as one with a mean cloud cover of less than 10 percent.

Code

##| out.height: 50%

    #prop clear days
    cloud_table %>%
      na.omit()%>%
      filter(month != 9)%>%
      mutate(binary_clear = dplyr::if_else(mean <= .1,true = 1,false = 0)) %>%
      group_by(id)%>%
      summarize(prop_clear = sum(binary_clear)/n(),
                clear_days = sum(binary_clear),
                total_days = n())%>%
      inner_join(x = boxes_sf)%>%
      ggplot(mapping = aes(fill = prop_clear))+
      geom_sf()+
      geom_sf(data = domain,inherit.aes = FALSE,fill=NA)+
      scale_fill_gradient(low = "white",high = "sky blue",limits=c(0,1))+
      geom_sf_text(aes(label = round(prop_clear,digits = 2)))+
      labs(fill = "Prop.\nClear",
           x=NULL,
           y=NULL)

Figure 5. Proportion of clear days (mean cloud cover < 10% ) for each flight box. Values in boxes are the proportion of clear days for that box

Diurnal patterns in wind speed

Code

era5_wind_table %>%
  inner_join(x = boxes_sf,
                 by = c("id"="ID"))%>%
  group_by(ordered_id, hour)%>%
  ggplot()+
  geom_point(mapping = aes(x=hour,y = wgt_90pct_wind_speed,color=target),alpha=0.5)+
  geom_smooth(mapping = aes(x=hour,y=wgt_90pct_wind_speed),method = "loess")+
  facet_wrap(~ordered_id,nrow = 2)+
  ylab("90th Percentile Wind Speed")+
  geom_hline(yintercept = 5, lty = 2)+
  scale_color_manual(values = c("aquatic" = "blue",
                                "terrestrial"="brown",
                                "Terrestrial (coincident aquatic)"="tan")) +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

Figure 6. Diurnal patterns in wind speed. Dots represent the 90th percentile wind speed in a flight box at a particular hour. Numbers above the plots refer to flight boxes. Lines shone are a loess fit.

Campaign Simulations

Code

  #code underlying simulations available in the file "R/flight_window.R"

      #readRDS("data/temp/sim_output_10pct.RDS")%>%
      readRDS("data/temp/sim_output_05pct.RDS")%>%
        #readRDS("data/temp/sim_output_01pct.RDS")%>%  
        #na.omit()%>%
        dplyr::filter(!is.na(box_id))%>%
        group_by(start_date)%>%
        summarise(sites_done = sum(!is.na(box_id)),
                  mean_cc = mean(na.omit(cloud_cover)),
                  days_taken = max(date)-min(start_date)+1)%>%
        ungroup()%>%
        mutate(day_of_year = yday(as_date(start_date)),
               year = year(as_date(start_date)))%>%
        group_by(day_of_year)%>%
        summarise(q0 = quantile(days_taken,probs=0),
                  q25 = quantile(days_taken,probs=0.05),
                  q50 = quantile(days_taken,probs=0.5),
                  q75 = quantile(days_taken,probs=0.95),
                  q1 = quantile(days_taken,probs=1)
        ) %>%
        mutate(date = as.Date(day_of_year,origin="2022-12-31"))%>%
        mutate(start_of_month = floor_date(date,unit = "month"),
               month_label = month(start_of_month,label = TRUE),
               julian_label = yday(start_of_month),
               day_of_month = mday(date),
               md_label = paste(month_label,"-",day_of_month,sep = ""))->test
      
      test %>%    
        ggplot()+
        geom_ribbon(aes(ymin=q0,ymax=q1, x = day_of_year),col="grey",alpha=0.2)+
        geom_ribbon(aes(ymin=q25,ymax=q75, x = day_of_year),col="grey",alpha=0.5)+
        geom_line(aes(x=day_of_year,y=q50))+
        geom_hline(yintercept = 37,lty=2)+
        scale_x_continuous(breaks = test$day_of_year,
                           labels = test$md_label)+ #note: it seems like it should be possible to inherit these
        ylab("Flight Days Needed")+
        xlab("Campaign Starting Day")+
        geom_text(label = "estimated total number of flight days",
                  y=37.3,
                  x=282)+
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) 

Simulations using cloud and wind data

In order to estimate how successful our campaign might be, we conducted a simulation using the MODIS aqua and terra cloud data.

1. Calculate mean cloud cover and median wind speed for each flight box for each day in the temporal window of interest (October-December, 2000 to present).

2. Rank flight boxes in descending order of priority. The aquatic sites receive the highest priority and are further ranked in descending order of median wind speed. The terrestrial sites are ranked in descending order of mean cloud cover. Thus, the highest priority is the windiest aquatic site and the lowest priority is the terrestrial site with the least cloud cover.

3. For each day in time time series:

1. If there are at least 60 days in which to sampling, continue. Else, skip.

2. Sample the highest ranking (i.e. hardest to sample) site (if any) that is below 5% cloud cover, has a median wind speed of less than 5 m/s(if it is an aquatic site), and hasn’t been sampled.

3. Repeat b until all sites have been sampled or no more time remains

The underlying code for the analysis is available at: https://github.com/BioSCape-io/terrestrial_sampling/blob/main/R/flight_sims_for_specific_days.R

Figure 7. Estimated length of time needed to finish campaign. The dark grey area represents the 50% CI, the light gray the 90% CI, and the solid line represents the median.

Figure 8. Fraction of campaign simulations in which a given box was flown.